Stochastic spreading processes on a network 
model based on regular graphs 



Sebastian V. Fallcrt^ and Sergei N. Taraskin^'^'* 

^ Department of Chemistry, University of Cambridge, Cambridge, UK 
^ St. Catharine's College, University of Cambridge, Cambridge, UK, 
sntlOOO@cam.ac.uk 

(* corresponding author) 



Abstract. The dynamic behaviour of stochastic spreading processes on 
a network model based on k-regular graphs is investigated. The contact 
process and the susceptible-infected-susceptible model for the spread of 
epidemics are considered as prototype stochastic spreading processes. 
We study these on a network consisting of a mixture of 2- and 3-fold 
coordinated randomly-connected nodes of concentration p and 1 — p, re- 
spectively, with p varying between and 1. Varying the parameter p from 
p = (3-regular graph of infinite dimension) to p = 1 (2-regular graph - 
ID chain) allows us to investigate their behaviour under such structural 
changes. Both processes are expected to exhibit mean-field features for 
p = and features typical of the directed percolation universality class 
for p — 1 . The analysis is undertaken by means of Monte Carlo sim- 
ulations and the application of mean-field theory. The quasi- stationary 
simulation method is used to obtain the phase diagram for the processes 
in this environment along with critical exponents. Predictions for critical 
exponents obtained from mean-field theory are found to agree with sim- 
ulation results over a large range of values for p up to a value of p = 0.95, 
where the system is found to sharply cross over to the one-dimensional 
case. Estimates of critical thresholds given by mean-field theory are found 
to underestimate the corresponding critical rates obtained numerically 
for all values of p. 
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1 Introduction 



The spread of epidemics poses a threat to biological populations as well as to 
computer networks and investigations into its dynamics and mechanisms are 
therefore of great current interest. One common class of epidemic models con- 
siders individuals to be in one of two possible states: susceptible (S) or infected 
(I). In this paper, we consider both the Contact Process (CP) [1] and the SIS 
model, two models of disease propagation via nearest neighbour contact, in which 



a disease is passed on to healthy nearest neighbours stochastically at a rate A 
specific to the model while infected sites spontaneously recover at rate e. 

These Markovian spreading processes have attracted wide attention in the 
past due to their applicability to phenomena as diverse as autocatalytic chemical 
reactions, spreading of rumours and transport in disordered media [2]. As the 
rates A and e are varied, an epidemic will be in one of two distinct states: an 
invasive regime (active state) in which it is present with a non-zero probability of 
ultimate survival and one in which this probability is zero thus leading to a state 
which allows no further evolution because the disease has died out (absorbing 
state). 

These two regimes are known to be connected by a continuous phase tran- 
sition thereby rendering them of conceptual interest for investigations into this 
kind of critical phenomenon of non-equilibrium statistical mechanics (see [3] 
for a review). The critical behaviour for these models in one-, two- and three- 
dimensional lattices has been investigated very accurately [4] and is found to be 
characteristic of the Directed Percolation (DP) universality class. From a range 
of studies, critical thresholds for the phase transition as well as critical exponents 
of predicted power-law scaling relations are known to high precision. 

With the growing interest in complex networks among the statistical physics 
community in recent years [5, 6], the question of the behaviour of dynamic pro- 
cesses on such topologically disordered structures has arisen [7, 8] . Particularly 
motivated by the fact that networked structures are ubiquitous in nature, the 
effects of these environments on, for example, the spread of a disease are of 
immediate interest. In a series of papers [9-13], the behaviour of the CP and 
the SIS model on a range of networks has been considered and even comparisons 
with data of computer virus outbreaks have been attempted [11]. As networks in 
general are infinite-dimensional objects, the dynamical mean-field (MF) approx- 
imation is expected to become exact in these cases in principle rendering many 
models tractable by analytical means. Both Monte Carlo (MC) simulations and 
the MF approximation have been used in previous investigations and produced 
such astonishing results as the absence of an epidemic threshold infection rate 
for infinite scale- free networks [11]. 

In this paper, we propose to investigate the behaviour of the CP and the 
SIS model as two paradigmatic stochastic spreading processes on networks of 
k-regular graph topology. The model network considered in this investigation 
consists of a mixture of 2- and 3-fold coordinated randomly-connected nodes of 
concentration p and 1 — p, respectively. Varying the parameter p from p = 
to p = 1 transforms the system from a 3-regular graph of infinite dimension to 
a 2-regular graph, i.e. a ID chain. While both the CP and the SIS model are 
expected to exhibit mean- field features for p — 0, the processes effectively take 
place in a one-dimensional environment for p = 1 which is a very well-studied 
regime of the DP universality class. It is our aim to investigate the behaviour of 
both the critical rates and accessible critical exponents for this crossover from 
an infinite- to a one-dimensional case thereby probing the validity of the MF 
approximation in this setting. The analysis is undertaken by means of Monte 



Carlo simulations using the quasi-stationary (QS) simulation method [14] and 
the application of mean-field theory [10]. 

This paper is structured as follows. Section 2 outlines the definitions and 
some properties of the processes considered. The MF approximation and the QS 
simulation method arc described in section 3. We present and discuss our results 
in section section 4 and summarise our findings in section 5. 

2 Background 

As outlined in the previous section, both the CP and the SIS model are simple 
toy models for the spread of an infectious disease by nearest-neighbour contact. 
In these models defined on a network, nodes represent susceptible or infected 
individuals surrounded by their neighbours connected via links along which the 
epidemic may spread. Proliferation of the disease to nearest neighbour sites hap- 
pens at a transmission rate A while recovery is spontaneous at rate e making the 
sequence of events an individual can cycle through Susceptible Infected 
Susceptible. 

The CP and the SIS model arc very similar, the difference being the exact 
mechanism of the spreading of infection. In the case of the CP, a site attempts 
to transmit its disease at rate A/fc to a randomly selected neighbour where k 
denotes the number of nearest neighbours. If the selected neighbour is already 
infected, proliferation fails. For the SIS model, transmission to any non-infected 
neighbour happens, in contrast, at rate A independent of the connectivity of the 
nodes. Thus, the spreading mechanism in the CP effectively compensates for 
the local connectivity present in the network through a suitable reduction of the 
spreading rate through a particular link. 

Once suitable initial states for all sites have been chosen, the above rules 
dynamically evolve the spread of a disease in the network. A typical initial 
condition is the state of a fully-infected system from which the system relaxes 
very quickly. For very long times, and formally as time t ^ oc and for an infinite 
number of sites in the network iV — > oo, the system is expected to be in one of 
two states: An active state in which there remains a finite density of infected 
sites or an absorbing state in which the disease has died out and that therefore 
admits no further time evolution. Depending on the value of the transmission and 
recovery rates A and e, ultimately the system will be in one of the two possible 
states. More precisely, there exists a continuous phase transition between these 
regimes as one fixes one of the rates and varies the other. This transition takes 
the system from a phase where the density of infected sites (order parameter) 
p is zero to one where it continuously grows from zero as the transmission rate 
(control parameter, assuming e fixed) is increased. 

Without loss of generality, one can perform a rescaling of time and set one 
of the two rates to unity for convenience. In the following, the recovery rate is 
assumed to be e = 1 and the critical point is therefore characterised by a critical 
transmission rate Ac alone. 



There exist a range of well-established scaling relations for various observ- 
ables in these models of which we present those relevant for this investigation. 
The density of infected sites in the thermodynamic limit as t ^ oo, the order 
parameter, is expected to scale as 

lim{p{t))=p^\X-Xef (1) 

thereby defining the order parameter critical exponent /3 where (...) denotes av- 
eraging over realisations of the process. Order parameter fluctuations are known 
to follow 

V = N(y^-f)^\X-X,r (2) 

Both the models under consideration are known to belong to the directed 
percolation (DP) universality class [2]. Accordingly, the critical exponents de- 
fined above arc those characteristic of this universality class. Above the upper 
critical dimension, d„ = 4, of these models, fluctuations are expected to be Gaus- 
sian and MF theory should be exact. Therefore, these processes taking place in 
infinite-dimensional networks are expected to exhibit exponents predicted by 
mean-field theory. 

3 Methods 

3.1 Mean-Field Approximation 

Both the CP and the SIS model can be described by the master equation which 
reflects the conservation of probability flow [3]. In the following, we will first 
outline the case of the SIS model and then consider the extension to the simpler 
case of the CP. 

In the dynamical MF approximation, which neglects density fluctuations and 
statistical correlations between the densities at different sites, and. for the mo- 
ment, disregarding the structure of the network completely, the master equation 
for the SIS model takes the form 

dtp{t) = -p{t) + A fc (1 - p{t)) p{t) (3) 

where p(t) denotes the density of infected sites at time t averaged over realisa- 
tions of the process which is identical to the probability of a site of the system 
to be infected at time t. This equation describes the rate of change of the aver- 
age density in the network which is equal to the flow of density into and out of 
any site with time and makes for the destruction and the creation terms above. 
The destriiction term due to the vanishing of infection at unit rate is propor- 
tional to the density p{t). The creation term is due to the possible infection by 
infected neighbouring sites in the case that the vertex under consideration is 
not infected. Accordingly, it is proportional to the probability that a site is not 
infected, (1 — p{t)), the probability that a neighbouring site is infected p{t), the 
local connectivity k and the spreading rate A. 



The master equation (3) can be extended in order to take into account the 
structure of the underlying network at the level of the node degree (connectivity) 
distribution (as developed by Pastor-Satorras and Vespignani [9]). It is clear 
that, unless one assumes a homogeneous network with (fc) ~ k for all k, the 
expression will dcciouplc into a set of equations for the densities of infected 
vertices characterised by a certain connectivity fc, We can write for each k, 

dtPk = -Pfe + A A; (1 - pk{t)) Ok{t) , (4) 

where 0k{t) is the probability that an edge emanating from a vertex of degree 
k is connected to an infected site. The infection term as described above now 
incorporates the probability that a site of degree k is connected to an infected 
vertex Ok{t). One can interpret Ok(t) as the mean density of neighbouring in- 
fected nodes and consequently k0k{t) as the mean number of infected nearest 
neighbours [9]. 

Networks which are Markovian are statistically described by their degree 
distribution P (k) and the conditional probabihty P{k'\k) that an edge of a 
node of degree k is connected to a vertex of degree k'. For such systems, we can 
write 

Ok{t) = Y, Pik'\k) Pk'it) (5) 

k' 

where the sum runs over all degrees k'. For uncorrelated networks, which we will 
exclusively consider in this investigation, this becomes [7] 



which docs not depend on k any longer. Substituting this expression into the 
rate equation (4) and imposing stationarity {dtPk = 0) one obtains 

Pk = = (7) 

where p^. and 0{\) are the time- independent values for the mean density at sites 
of degree k and the mean density of infected neighbouring sites. Multiplying by 
^-^1^ and summing over k yields 

1 l^P(fc)_fc^^ 

where f{x) is a monotonically decreasing fmiction of x. This equation only has a 
(unique) solution for 6* (A) different from zero for A > Ac where Ac is the threshold 
value for the transmission rate that makes O smallest. As by definition the mean 
density of infected nearest neighbours > in the active regime and f{x) is a 
monotonically decreasing function we have (effectively setting = 



for the critical threshold [15]. 

In order to obtain an expression for the order paramter we start by combining 
equations (6) and (7) and arrive at a self- consistency equation for (equivalent 
to equation (8)) 

y^m^J^ (10) 

^ {k) i + k\e ^ ' 

that can in principle be solved for which in turn allows one to obtain the order 
parameter in the MF approximation from 



Critical exponents can be extracted from MF theory by considering the lead- 
ing behaviour of the relevant expressions. For example, combining the last ex- 
pression Eq. (11) and Eq. (7) and expanding in A — in analogy to the scaling 

form p ~ (A — Ac)'', one obtains the MF value for the order parameter exponent 
(5=1. Similarly, one obtains 7 = for the corresponding fluctuations. 

In the case of the CP where the effective spreading rate is inversely propor- 
tional to the number of links connected to an infected site, Eq. (5) has to be 
modified and reads 

k' 

which for uncorrelated networks leads to 0"''^{t) = p{t)/{k) where p is the average 

density of infected sites averaged over degrees. Following the procedure as for 
the SIS model, the critical threshold rate is found to be degree independent and 
given by [12] 

A?P = 1 (13) 
while the critical exponents are identical to those for the SIS model. 



3.2 Monte Carlo Simulation: Quasi-Stationary Simulation 

Both processes under consideration can be simulated effectively via time-dependent 
Monte Carlo simulations. Once initial conditions have been chosen, the system 
is evolved according to the appropriate rules with a simulation that selects pos- 
sible events according to their prescribed rates and ensures that they happen in 
exponentially-distributed time intervals. 

In contrast to lattices, networks are characterised by the existence of long 
range links. This implies that a dynamical process will very strongly feel the size 
of the system in a finite representation of the network leading to much stronger 
finite-size effects than experienced in lattices. While the most precise method of 
determining the critical point in lattices is spreading from a single seed while 



ensuring that the infection never reaches the boundary of the system, this is 
virtually impossible in networks. Thus, finite-size eS'ects have to be systemat- 
ically exploited in order to make a prediction about the infinite system using 
networks of different sizes. This is possible using the finite-size scaling hypothe- 
sis which predicts that values of obscrvablcs in systems of size L are controlled 
by the ratio L/^±, where is the spatial correlation length. In order to use this 
scaling behaviour, one requires a stationary average value of these observables 
and analyses their values for different system sizes. Problematically, due to the 
existence of the absorbing state no such true stationary state exists in a finite 
system. Fortunately, the processes under consideration evolve such that some 
observables attain quasi-stationary (QS) values. 

Most notably, the density of infected sites averaged over siurviving realisa- 
tions of the process {p{t)) exhibits this behaviour after an initial transient start- 
ing from the initial state of a fully- infected system [2]. This quasi-stationary 
density p is expected to scale systematically in accordance with the finite-size 
scaling hypothesis. According to recent investigations into the finite-size scaling 
behaviour above the upper critical dimension [16, 13] the prediction is 

p - L^f"/^ G (^L^/^X - Ae)) = 7V^/2 Q (7vi/2(A - A^)) (14) 

where L"^ ~ N, the number of sites in the system and G{x) an appropriate scaling 
function. A similar expression with /? replaced by 7 is valid for the associated 
fluctuations. 

In principle the QS state can be investigated via a conventional simulation 
in which the system is stochastically evolved in time from a fully-infected initial 
condition. The density of infected sites conditioned on survival can then be 
analysed and a temporal average over the duration of the QS state is an estimator 
for p. This method is however plagued by a range of problems [17] which led de 
Oliveira and Dickman to propose a simulation method which samples the QS 
state directly [18]. 

In this QS simulation method, the absorbing state is eliminated and its prob- 
ability weight is redistributed over the active states according to the history of 
the process. It can then be shown that the true stationary state of the resulting 
modified process corresponds to the quasi-stationary state of the original one. 
This method is ideally suited for our study as a single realisation of the network 
is investigated in one QS simulation run enabling us to analyse samplc-to-sample 
fluctuations between realisations and find a critical point by use of the scaling 
relation Eq. (14). 

4 Results and Discussion 

4.1 MF Solution 

The mean-field theory for the CP and the SIS model on networks can be applied 
to our network with two types of nodes of connectivity ki and ^2 present with 
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Fig. 1. The mean-field average density of infected sites for tlie SIS model on a binary 
random network with fci = 2, ^2 = 3 for (from top to bottom) p=0, 0.5, 1. 



probabilities P{ki) and P(k2) respectively. For the SIS mode, the self-consistency 
equation, Eq. 10, for the stationary value of 0, the average density of infected 
nearest neighbours, takes the form 



A 



1 + kixe 1 + fczAe 



which can be solved for O giving 



e 



1 

2A 




ki +k2 
kik2 



{k^)X - (k) 
kik2{k) 



(15) 



(16) 



where is only defined for transmission rates A > Ac = as explained 
above. Using the definition of the average stationary density p — J2k ^i^)Pk 
and finally substituting P{ki) — p and P{k2) = 1 — p, the order parameter in 
the MF approximation is given by 



\0 



pki 



I (1-P)fc2 

l + kiX0 l + k2X0 



(17) 



The critical exponents given by the leading order contributions of the relevant 
expressions are found to be the standard MF exponents [12] as outlined above. 

This solution is plotted for the special case ki = 2 and ^2 = 3 with p varying 
from to 1 in Fig. 1. As expected, MF theory qualitatively reproduces the 



features of the phase transition: There exists a critical rate Ac below which the 
stationary density is zero while it grows continuously for values above. For the 
CP, an analogous analysis yields a similar solution. 

No direct comparison with numerical predictions for the order parameter 
is shown in Fig. 1 because of severe finite-size effects for networks which shift 
the simulation curve well above its true asymptotic position. We will however 
compare the critical threshold rate as well as some critical exponents in the next 
section. 

4.2 Simulation Results 
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Fig. 2. The critical threshold Ac obtained from MP theory (lines) and MC simulation 
(circles and squares) for the CP (circles and solid line) and the SIS model (squares and 
dashed line) for various values of p. Inset shows the critical exponents j3 (circles) and 
7 (squares) as a function of p. 

The critical thresholds Ac for a range of values of p for both the CP and the 
SIS model were obtained via the QS simulation method outlined above. We used 
networks of sizes ranging from = 256 — 32768 in QS simulation runs up to 10^ 
time steps averaging over no less than 100 and up to a maximum of 1000 network 
realisations (the latter were required to minimise errors in light of strong sample- 
to-sample fluctuations for large values of p) . The resulting critical thresholds are 
shown in Fig. 2 along with the MF predictions obtained from Eqs. (9) and (13). 
As expected from the definition of the two models, the CP threshold exceeds the 
one of the SIS model for a particular value of p which can be attributed to the 



reduction of the effective transmission rate by the local coordination number as 
in Eq. (12). Also, in the two cases of homogeneous connectivity, p — and p = 1, 
the thresholds for the two models are expected to be simply related by a factor 
of 3 and 2 (equal to the connectivity), respectively, as can be seen from Fig. 2. 
Note that the critical threshold for the SIS model in the quasi one-dimensional 
case (p = 1), Ac = 1.65, is almost identical to the threshold for the CP on the 
3-regular random network {p = 0) Ac = 1.63. 

Turning to the MF predictions in comparison to the MC results, one notes 
that they underestimate the true critical thresholds for all values of p and for 
both models. The difference between the MF approximation and the simulation 
results is more pronounced for the CP as compared to the SIS model. 
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Fig. 3. The density of infection p (upper panel) and the corresponding fluctuations 
V = N (^p^ — (lower panel) for various network sizes for the case p — 0.5 in the SIS 
model. Solid lines are best-fit regression lines to the scaling forms defined in Eq. (14). 

The exponents (3 and 7 were obtained by fitting data for the density of 
infected sites p and the corresponding fluctuations V — N — to the 

finite-size scaling form of Eq. (14). A typical set of data points is shown in Fig. 3 
for the case p = 0.5 for the SIS model along with best-fit regression lines. Both 
quantities as a function of network size N show power-law behaviour with the 
expected MF exponents (3/2 = 0.5 and 7 = indicating the validity of the MF 
approximation for this case. These values of critical exponents are plotted in the 
inset of Fig. 2. As p is further increased, strong sample-to-sample fluctuations 
for values beyond p = 0.95 render a precise analysis very complicated. For p = 1, 



the well-established ID finite-size scaling exponents are recovered {I3/v_l = 0.253 
for p{N) and = 0.498 for V"(A'') [2]) as can be seen from the figure. In 

the transition region, error bars for exponents are large and our results give 
slight preference to the scenario of a discontinuous change in exponent values. 
However, we feel that a very rapid yet continuous change of exponents towards 
the ID values cannot be excluded. 

5 Conclusion 

Wc have investigated the CP and the SIS model, two paradigmatic stochastic 
spreading processes, in a network model which interpolates smoothly between 
an infinite-dimensional 3-regular random network and a linear chain through 
the variation of a single parameter p. The MF approximation yields a prediction 
for the critical threshold rate of an epidemic outbreak and critical exponents 
associated with the corresponding absorbing state phase transition. For no value 
of p docs MF theory predict the true critical threshold as calculated from MC 
simulations. The predictions for critical exponents agree perfectly with simula- 
tions for a very wide range of p up to p = 0.95. Beyond this point, the analysis 
is complicated by strong samplc-to-sample fluctuations. For p = 1 one recov- 
ers the established exponents for the one-dimensional case indicating a sudden 
crossover. While not being able to investigate the nature of this transition pre- 
cisly, our simulations favour the scenario of a discontinuous change in the scaling 
exponents which reflects the abrupt change of the dimensionality of the network. 
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